subroutine init_ch()
  use mod_para
  use mod_workspace
  use mod_material
  use mod_grid
  use mod_migrate
  implicit none
  integer::i,j,k,sourcex,sourcey
  real(8)::ephifactor,r,tstar,zz,ephi,ephimax
  

  real(8)::t1,t2,t3
  real(8)::a1x,a1y,a1z,a2x,a2y,a2z,a3x,a3y,a3z
  real(8)::ca1x,ca1y,ca1z,ca2x,ca2y,ca2z,ca3x,ca3y,ca3z,posx,posy,posz
  


  sourcex=(nx1+nx2)/2
  sourcey=(ny1+ny2)/2
  posx=(x(sourcex,sourcey,nz1)+x(sourcex+1,sourcey,nz1))/2D0
  posy=(y(sourcex,sourcey,nz1)+y(sourcex,sourcey+1,nz1))/2D0
  posz=z(sourcex,sourcey,nz1)

  tstar=1.13D0*delta**2
  ephimax=-1D100

  do k=nz1,nz1+1       
     ephifactor=1/(2.D0*pi*sigma(sourcex,sourcey,nz1))
     do j=ny1,ny2
        do i=nx1,nx2
           zz=z(sourcex,sourcey,k)-z(sourcex,sourcey,nz1)
           r=sqrt( (x(i,j,k)- posx )**2+ ( y(i,j,k)-posy )**2)+epsilon(1D0)      
           call calcephi(r,zz,tstar,ephi)
           ephi=ephi*ephifactor
		   ephimax=max(ephi,ephimax)
           ey(i,j,k)=ephi*(x(i,j,k)-posx)/r
           ex(i,j,k)=-ephi*(y(i,j,k)-posy)/r      
        enddo
        write(*,*)  'j=', j,' k=',k,'ephimax=',ephimax
		ephimax=-1D100
     enddo
	 arrayview=ex(:,:,k)
  enddo



  do k=nz1,nz2
     do j=ny1,ny2
        do i=nx1,nx2
           a1x=x(i+1,j,k)-x(i,j,k)
           a1y=y(i+1,j,k)-y(i,j,k)
           a1z=z(i+1,j,k)-z(i,j,k)
           a2x=x(i,j+1,k)-x(i,j,k)
           a2y=y(i,j+1,k)-y(i,j,k)
           a2z=z(i,j+1,k)-z(i,j,k)
           a3x=x(i,j,k+1)-x(i,j,k)
           a3y=y(i,j,k+1)-y(i,j,k)
           a3z=z(i,j,k+1)-z(i,j,k)
 

           ! compute covariant ex from ordinary ex
           t1=( a1x*ex(i,j,k)+a1y*ey(i,j,k)+a1z*ez(i,j,k)  )
           t2=( a2x*ex(i,j,k)+a2y*ey(i,j,k)+a2z*ez(i,j,k)  )
		   t3=( a3x*ex(i,j,k)+a3y*ey(i,j,k)+a3z*ez(i,j,k)  )
           ex(i,j,k)=t1
           ey(i,j,k)=t2
		   ez(i,j,k)=t3
        enddo
     enddo
	 arrayview=ex(:,:,k)
	 arrayview=ez(:,:,k)
  enddo

  
  do k=nz1,nz2
     do j=ny1,ny2
        do i=nx1,nx2
           chx(i,j,k)=(ey(i,j,k+1)-ey(i,j,k)+ez(i,j,k)-ez(i,j+1,k))/(mu*jac(i,j,k))
           chy(i,j,k)=(ex(i,j,k)-ex(i,j,k+1)-ez(i,j,k)+ez(i+1,j,k))/(mu*jac(i,j,k))
           chz(i,j,k)=(ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k))/(mu*jac(i,j,k))
        enddo
     enddo
	 arrayview=chx(:,:,k)
	 arrayview=chz(:,:,k)
  enddo

  chz(:,:,nz2+1)=0

  do k=nz2,nz1,-1
     do j=ny1,ny2
        do i=nx1,nx2
           ! chz(i,j,k)=-chx(i,j,k)+chx(i+1,j,k)-chy(i,j,k)+chy(i,j+1,k)+chz(i,j,k+1)
           chz(i,j,k)=( &
                -chx(i,j,k)*jac(i,j,k)+chx(i+1,j,k)*jac(i+1,j,k) &
                -chy(i,j,k)*jac(i,j,k)+chy(i,j+1,k)*jac(i,j+1,k) &
                +chz(i,j,k+1)*jac(i,j,k+1) &
                )/jac(i,j,k)
        enddo
     enddo
	 arrayview=chz(:,:,k)
  enddo

end subroutine init_ch



